
options(warn=1)
args = commandArgs(trailingOnly=TRUE)
if (length(args)!=3) {
  stop("Missing arguments: Provide DD threshold and path to split 0 and path to split 1")
}

kDDThreshold <- as.numeric(args[1])
kPathSplit0 <- as.character(args[2])
kPathSplit1 <- as.character(args[3])

set.seed(7)

### each dataset consists of two rows (one for X and one for 1-X); here, we only consider 1 row for each dataset (results are similar for the other row or both rows)
evenOdd <- "odd"
#evenOdd <- "even"
#evenOdd <- "all"  



totalDatasets <- 459 #
splitMask <- numeric(totalDatasets) + 1 


# skip=1  The first row is 'summary'
res=read.csv(file=kPathSplit1, skip=1, header=TRUE, sep=",")
res2=read.csv(file=kPathSplit0, skip=1, header=TRUE, sep=",")

if (evenOdd == "odd") {
  res=res[seq(1,nrow(res),2),]#selects odd indexes for more independent capture of the proposed bound.
  res2=res2[seq(1,nrow(res2),2),]
} else if (evenOdd == "even") {
  res=res[seq(2,nrow(res),2),]#selects even indexes for more independent capture of the proposed bound.
  res2=res2[seq(2,nrow(res2),2),]
}
stopifnot(length(c(res[,1], res2[,1]))==totalDatasets)

uid=c(res[,1], res2[,1])[which(splitMask==1)] 
b=c(res[,2], res2[,2])[which(splitMask==1)] 
bl1=c(res[,3], res2[,3])[which(splitMask==1)]  
bu1=c(res[,4], res2[,4])[which(splitMask==1)]  
cl=c(res[,5], res2[,5])[which(splitMask==1)] 
cu=c(res[,6], res2[,6])[which(splitMask==1)] 
dl=c(res[,7], res2[,7])[which(splitMask==1)] 
du=c(res[,8], res2[,8])[which(splitMask==1)] 
t3=c(res[,9], res2[,9])[which(splitMask==1)] 
sst=c(res[,10], res2[,10])[which(splitMask==1)] 
ssreslo=c(res[,11], res2[,11])[which(splitMask==1)] 
ssres2=c(res[,12], res2[,12])[which(splitMask==1)] 
p=c(res[,13], res2[,13])[which(splitMask==1)] 
l=c(res[,17], res2[,17])[which(splitMask==1)] 
u=c(res[,18], res2[,18])[which(splitMask==1)] 
ssreslob=c(res[,23], res2[,23])[which(splitMask==1)] 
ssres1b=c(res[,24], res2[,24])[which(splitMask==1)] 
#t1b=res[,25]
tb2=c(res[,31], res2[,31])[which(splitMask==1)] 
tw2=c(res[,32], res2[,32])[which(splitMask==1)] 
tb1=c(res[,33], res2[,33])[which(splitMask==1)] 
tw1=c(res[,34], res2[,34])[which(splitMask==1)] 

datasetLabel <- c(as.character(res[,36]), as.character(res2[,36]))[which(splitMask==1)] 


###########
bud=pmax(pmin(bu1,du),dl)
bld=pmin(pmax(bl1,dl),du)
sl=bl1-cl
su=cu-bu1
###############################

xs=c(0.00,0.25,0.50,0.75, 1.00, 1.25, 1.50, 1.75, 2.00)
ms <- numeric(length(xs))
my <- numeric(length(xs))
mw <- numeric(length(xs))
conf <- numeric(length(xs))


for (j in 1:(length(xs) )) { 

  ex=xs[j]
  
  bu=bu1+ex*su
  bl=bl1-ex*sl
  bbl = pmax(bl,dl) 
  bbu = pmin(bu,du) 
  
  #
  y=1*(b<  bbu )*1*(b> bbl )
  wr=(bbu-bbl)/(du-dl)
  
  sel=1-1*(pmin(bu1,du)<pmax(bl1,dl))
  
  distanceSel <- numeric(length(sel))

  for (i in 1:length(distanceSel)) {
    if ((du[i]-dl[i] < kDDThreshold)) {
      distanceSel[i] <- 1
    }
  }
  sel=sel*distanceSel
  
  ms[j]=mean(sel)
  my[j]=mean(y*sel)/mean(sel)
  mw[j]=mean(wr*sel)/mean(sel)
  conf[j]=1-pnorm(-ex)

  mind=pmin(abs(b-bbl)[which((1-y)*sel==1)],abs(bbu-b)[which((1-y)*sel==1)])
  md=mean(mind)

  x=ex
}
print(sprintf("Number of datasets considered: %d", length(uid)))
print(sprintf("Selected percent: %f", ms[3]))

cat("$x$ & $\\Phi(x)$ & p($B \\in CI_x | selected)$ & $E[WR_x | selected]$ \\\\ \n")
for (printIndex in 1:length(xs)) {
  cat(sprintf("%s & %s & %s & %s\\\\ \n", format(round(xs[printIndex], 4), nsmall=2), format(round(conf[printIndex], 4), nsmall=4), format(round(my[printIndex], 4), nsmall=4), format(round(mw[printIndex], 4), nsmall=4)))
} 
